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^^ ' We study the eflfect of shear flow on homogeneous crystal nucleation, using Brownian Dynamics 

^N ^ simulations in combination with an umbrella sampling like technique. The symmetry breaking due 

f~| ■ to shear results in anisotropic radial distribution functions. The homogeneous shear rate suppresses 

'—{ ' crystal nucleation and leads to an increase of the size of the critical nucleus. These observations can 

y~~i , be described by a simple, phenomenological extension of classical nucleation theory. In addition, we 

I . find that nuclei have a preferential orientation with respect to the direction of shear. On average 

I ' the longest dimension of a nucleus is along the vorticity direction, while the shortest dimension is 

I preferably perpendicular to that and slightly tilted with respect to the gradient direction. 

O ' I. INTRODUCTION 

C^ , It is well known that fluids can be cooled to temperatures below coexistence temperature of the liquid and solid 

G ' phase. Density fluctuations in this super-cooled regime lead to the continuous formation of small crystalline nuclei. 

I , Even though the crystalline phase is the thermodynamical stable phase, a super-cooled system can remain for longer 

' O ■ times liquid-like. The reason lies in the fact that there is a counter acting force to crystallisation which finds its 

origin in the interface that needs to be formed between small crystallites and the surrounding liquid. This results in 

a nucleation barrier which needs to be overcome in order to solidify the system. This intuitive picture forms the basis 

of what at present is known as classical nucleation theory [l[ y] • 

Understanding the underlying principles of nucleation and the growth to macroscopic sizes finds its applications in 
^^ ' many areas, e.g. protein crystallisation and metallurgy |^y,|5|. The first experimental studies of crystal nucleation 
J^ are light scattering experiments 6, 7j,||8j in hard sphere colloids, which form a well understood model system |9|. Crys- 

r — ' tallisation rates were measured by Gasser et al. using confocal microscopy in weakly charged colloids |l(l |. Recently 
ps) [ computer simulations have been used to explore nucleation phenomena and by using special simulation techniques it 
VO ' was possible to obtain quantitative results on the absolute rate of crystal nucleation in colloidal suspensions [llL lia | 
GD I and compare them with experiments. 

^ . In the present work we focus on nucleation outside equilibrium, i.e. in the presence of a homogeneous shear 

field. Experiments on crystallisation under shear suggest two different scenarios. In some cases one observes that 
crystallisation is enhanced by the presence of shear |lv3 , J[4 , J^, jl6,] . This can, at least qualitatively, be understood, 
by the idea that shear might facilitate the formation of layers in the system, making it easier to crystallise. Other 
experiments report that crystallisation is suppressed [13, ■ In the later case one could imagine that the shear will 
i-rt ' de-attach particles from newly born crystallites. Both effects are also found in heterogenous nucleation under shear 
Pj ' in single or double wall confinement. For low and moderate shear rates two different types of layering sliding are 
O 1 observed, while for high shear rates the layers are found to destabilise [l9j . 

^ ' A priory it is not clear which of the both scenarios will dominate and how it depends on the external conditions. In 

J> simulations by Butler and Harrowell '2(J| for the crystallisation kinetics of particles with Yukawa interaction subjected 

to shear rates, it was observed that crystallisation is suppressed by shear. In the present work we will focus on 

nucleation in its bare bones and extract information on the properties of nuclei that are formed in a super-cooled 

j^ ■ melt. A preliminary account of our simulation results was already published |2l|. We include here further data for 

pair correlations, the critical nucleus structure and orientation, and describe the simulation technique in more detail. 

The remainder of this paper is organised as follows. In Sec^Jwe introduce the model and explain the combination 

of different techniques we used for our simulations. The simulation results for the pair-correlation functions and the 

size distribution function of the nuclei are shown in Sec. lIIII In Sec. II VI we make a simple, phenomenological extension 

of classical nucleation theory and compare this with the results obtained from the simulations. A more detailed 

analysis of the shape and relative orientation with respect to the shear direction is performed in Sec. We finish in 

Sec lVII with a summary and discussion of our main results. 
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II. THE MODEL 

In order to simulate homogeneous crystal nucleation under shear, we have chosen to consider the simple model 
of charge-stabilised colloidal suspensions subjected to a linear shear flow. Due to screening of the micro-ions, these 
particles mutually interact effectively via a repulsive Yukawa potential Q 

p-nr 

yW = e^, (f) 

where r is the distance between the centres of the colloidal particles. The magnitude of screening is expressed by the 
inverse screening length n and the strength of the interaction is denoted by e. 

In addition to the normal inter-particle interaction, we also want to apply shear to the system. This means 
that we are dealing with a non-equilibrium situation, hence we need to make use of dynamical simulations rather 
than simulations of the Monte Carlo type. For low volume fractions of charged suspensions one can neglect the 
hydrodynamic interactions between the colloids, which enables us to apply Brownian Dynamics |23 . l2a| . In this 
over-damped motion, it is assumed that particles move at constant velocity proportional to the force due to other 
particles, and that the solvent exerts a random force with a Gaussian distribution. The steady shear rate 7 imposed 
on the system, is assumed to generate a linear velocity gradient in the solvent, which in turn results in a linear force 
field for the colloids. The combination of Brownian Dynamics and shear results in the following equations of motion 

n{t -f St) = n{t) + St^ + 6r^ + St^y,{t)x . (2) 

Here ri(i) = (a;i(t), j/i(t), Zi(i)) is the position of the ith colloidal particle at time t. In a small time interval 6t 
this particle moves under influence of the sum of the conservative forces fi{t) arising from the pair interaction |^ of 
particle i with the neighbouring particles. The friction constant ^ with the solvent is related to the short-time diffusion 
constant D by ^ = /cbT/Z). The stochastic displacements are independently drawn from a Gaussian distribution with 
zero mean and variance {{Sr^^)"^) = 2D5t, where a stands for one of the Cartesian components. The last term in 
Eq. ^ represents the applied shear in the x-direction, and imposes an explicit linear flow field. 

All simulations reported here are done for a system of 3375 colloidal particles and the interaction strength is fixed 
to e = 1.48 X W^ksT. For practical purposes the interaction potential is truncated at the cut-off distance re = 10/k, 
by shifting the potential to zero. The inter-particle forces are not effected by this procedure. For simplicity we use 
a cubic simulation box to which we apply the shifted periodic boundary conditions. These were introduced by Lees 
and Edwards [23| and are required to deal with the presence of shear. 

The simple combination of Brownian Dynamics with a linear imposed velocity field and Lees-Edwards boundary 
conditions, forms the heart of the simulations. This ensures that the starting configurations follow the proper equations 
of motions. To facilitate the connection with Gibbs free energies, however, it is desirable to do the simulations 
at constant osmotic pressure. Hereto we have added volume moves as normally used in equilibrium Monte Carlo 
simulations of the isobaric ensemble. In practice this means that after a large number of Brownian Dynamics steps, 
the simulation box expands or shrinks isotropically. The difference in potential energy due to this volume move, is 
used in order to accept or reject this move. The value of the applied pressure is chosen such that the bulk behaviour of 
the system is in the super-cooled regime. Therefore the prepared system will remain liquid, even though the crystalline 
phase is more stable. Due to fluctuations inside such a liquid, small nuclei are continuously formed and dissolved. 
Only once a sufficiently large nucleus is formed, the critical nucleus size, the growth of the nucleus will dominate and 
can lead to a full crystallisation of the system. 

The quantity of interest of our simulations is the cluster size distribution function P{n), which describes the 
probability of finding a nucleus of n particles. This is a well defined object, even in the non-equilibrium situation 
of the presence of shear, provided that after a sufficient long simulation time the system will reach a steady state 
situation. To this end we need to analyse the configurations obtained by the simulations and detect all the nuclei in 
the liquid. This can be done with the aid of bond-orientational order parameters introduced by Steinhardt et al. |25| 
and applied to study crystal nucleation by Frenkel and coworkers ,2£, .22^ ^] . All particles within a short distance of 
the particle of interest are determined, and the vectors connecting them are resolved in polar angles. These can be 
used to determine rotationally invariant order parameters, which characterise the neighbourhood of a particle, and 
allow us to distinguish a local fluid environment from a crystalline environment. It is even possible to distinguish 
between different crystalline structures. Solid particles that are in each others neighbourhood belong to the same 
cluster. In this manner we obtain the total number of clusters in the system and each of their sizes. 

By performing long simulations, in which between 10^ and 10^ Brownian Dynamics steps are made, we can now 
measure the probability size distribution function P{n) in a simple way. After 20 Brownian Dynamics steps and 
possibly a volume move (on average one per hundred BD steps), we analyse the configuration in determining the size 



and the number of clusters present in the hquid. These are used to make a histogram of the possible cluster sizes, 
and the frequency is proportional to the probability of observing such a cluster size. Note that also in the presence of 
shear, and therefore in non-equilibrium, this procedure can be followed. 

Nucleation, however, is a rare process and unless the system has a large nucleation rate :.29], no sufficient statistics 
can be obtained by merely the technique described above. The solution to this problem is the usage of umbrella 
sampling |23 • This technique stems from equilibrium Monte Carlo simulations and introduces an additional potential 
to the system. Such a bias should depend on a relevant reaction coordinate of the system under consideration and 
is able to restrain the system in meta-stable situation. The traditional interpretation is that the bias is used as an 
external potential, leading to a biased ensemble. It therefore might seem somewhat surprising that the technique can 
be used even in non-equilibrium situations. The reason for this apparent ambiguity, is that umbrella sampling can be 
considered as a plain mathematical trick. 

In order to measure the steady state probability distribution P{n), we need to sample the relative probabilities of 
P{n~\- 1)/P(n), which is achieved by growing paths by Brownian dynamics simulation from one cluster size to another. 
In principle this could be done by starting with a chosen cluster-size n and count the number of times the cluster 
grows or shrinks with one particle or remains at the same size in a small interval of time, and average afterwards over 
many different realizations of the initial configuration. The number of counts of a given size n is then proportional 
to P{n). The obvious problem is that below the critical nucleus size, the cluster will mainly shrink, and beyond 
the critical size it will mainly grow. This would be alright, but causes a different problem, i.e. how to obtain many 
configurations of a given size which have a very low probability. 

This problem is circumvented by using the bias. For each simulation we use a fixed biasing function based on some 
preferred cluster-size no, which has the simple form 

Vbias{n) = -a[n - tiq)^, (3) 

with a is a measure for the strength of restrainment. We start with a given configuration and follow the normal 
Brownian Dynamics as described above. After a fixed time interval we measure the size n of the largest cluster. We 
now continue with either the final configuration or the previous starting configuration as the new initial configuration. 
Only this is determined by the bias we introduce, it therefore has no effect on the dynamics of the system and is just 
an aid to generate new initial configurations. By choosing the bias appropriate to a harmonic-like form, we ensure 
that the probability of growing a bigger or shrinking to a smaller cluster with respect to no is approximately the 
same. The count of the number of clusters of a given size, is now proportional to the ratio of P(n) and PuasifT') of 
the known bias we impose. By performing simulations for different preferred sizes ng, we obtain successive parts of 
the function P{n) which can be combined to form the complete function. 

III. SIMULATION RESULTS 

The pressures at which simulations are performed are such that the equilibrium bulk phase without shear is a 
face-centered-cubic crystalline phase |12t i31j . Hence the liquid state in which we prepare the system is only meta- 
stable. Nevertheless, the pressure is still small enough such that spontaneous crystallisation is not observed within 
the duration of the simulation and even in the absence of the bias the system remains liquid and only forms nuclei 
smaller than the critical nucleus size. 

One of the consequences of the presence of shear is that this introduces a symmetry breaking in the fluid. This can 
be seen easily in the 2-particle distribution function g^'^'{fi,f2)- Whereas in bulk-equilibrium, this function is only 
dependent on the absolute distance |ri — r2|, in the presence of shear the function will also depend on the relative 
direction. In general the function can be expanded in a series of modified spherical harmonics C;^m(0, (p) with distant 
dependent coefficients gi^m{r) 

°° ' 21 + 1 

where 9 is the angl e with resp ect to the vorticity direction, i.e. cos6' = f • z, and the angle with the shear direction, 
i.e. cos (p = f ■ xj \J\ — (f • zY (See Fig.^l. Note that without the symmetry breaking due to shear, only the isotropic 
contribution remains, i.e. the usual radial distribution function g(r) = go.o{r). 

The distance dependent coefficients follow directly from the orthogonality relation of the spherical harmonics 

gi,rnir) = J dcos9 J d.p gir,9,(P)Cl^i9,^) (5) 
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FIG. 1: Schematic drawing of the different directions in the problem. The reman numbers label the different quadrants of 
space. 

and can easily be obtained from simulations. Due to symmetries present in the system and the suitable chosen 
reference frame, only functions with both indices I and m even have non-zero contributions. It turns out that the 
first order correction to the radial distribution function is dominated by the imaginary part 3(g2,2(^)) of the function 

52,2(0- 

In Fig.[5]we show the ratio of 3(g2.2(^)) and the function gQ.o{r) as obtained from simulations without bias and for 
different shear-rates. The corresponding function, which was statistically averaged, is S(C2,2) = — |:-\/3sin^ 6'sin2(/). 
This quantity is able to distinguish whether a particle is preferably found in the first (third) or in the second (fourth) 
quadrant of space provided another particle is in the origin (See Fig.^. The signal is phase-shifted with respect to 
the oscillations in the radial distribution function, such that particle distances in a peak of a distribution function 
and located in the first quadrant are typically closer than those particles found in the second quadrant. It can also 
be observed that the amplitude grows linearly with the shear rate. The distribution functions 50.0 (^) are, within the 
accuracy of the simulations, identical for these moderate shear rates. 



0.02 



0.01 



o) 0.00 



S -0.01 



-0.02 



-0.03. 



1 1 1 1 1 1 i 




/Tf Tf ^^^^ 


- 


— t/k'd = (g(r)) 




'- — JIkD = 0.0008 


- 


1 — y/K'D = 0.0016 




' — y/K"D = 0.0032 
I.I.I, 



10 



20 
rK 



30 



40 



FIG. 2: The imaginary part of the function 32,2(7') as function of the distance for fixed pressure PP/ k. — 0.240 and different 
shear rates. For comparison the radial distribution function g{r) is included. 



In order to obtain the probability function P{n) of observing a cluster of n particles, we perform simulations for 
fixed pressure and shear-rate. The values of the preferred size no in the bias Q is chosen in multiples of 10 and the 
strength a = 0.15. This ensures that the window of observed cluster sizes are sufficiently wide to produce overlap 
between neighbouring windows. When the bias is included, we do not only find a single big nucleus, but also small 
nuclei can be observed. It suffices, however, only to count the largest cluster, as it are those fiuctuations in size that 
determine the probability function P{n). The overlap between neighbouring windows is used to match successive 
parts of P{n). The only exception is the first part of the function, where small nuclei are important. For that region 



an unbiased simulation is used in which all clusters in the system are used. For a more detailed discussion on how 
the curve is obtained we refer the reader to |28| . 

The results of the simulations for a single pressure and several shear-rates are shown in Fig. 13 It shows the negative 
logarithm of the relative probability function P{n) normalised with the probability P{1) of finding a cluster of only 
a single particle. For increasing values of the shear rate, the maximum in the curves shifts towards larger values and 
larger cluster sizes. 
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FIG. 3: Negative logarithm of the probability P{n) of finding a cluster of n solid-like particles normalised by P(l) for pressure 
/3P/k3 = 0.24 and different applied shear rates, from bottom to top -y/iK^D) = 0, 0.8 x 10"^, 1.6 x 10"^, and 3.2 x 10"^. The 
insets show typical snapshots of critical nuclei for the largest shear rate and the zero shear case. 



IV. EXTENSION OF CLASSICAL NUCLEATION THEORY 

Classical nucleation theory links the probability of finding a cluster of given size to the Gibbs free energy and allows 
to interpret a probability function as shown in Fig. O for the zero-shear case in terms of nucleation barrier [23 . It 
splits the cost in Gibbs free energy AG to create a spherical nucleus in a liquid in two contributions 



AG: 



-irR'^PslAfil+ATrR^-fSL- 



(6) 



On the one hand there is a gain in energy proportional to the volume of the nucleus with radius R due to the difference 
in chemical potential A/i between the solid with density ps and the liquid phase. On the other hand we have a loss 
in energy, since an interface between the solid nucleus and surrounding liquid needs to be formed, described by 75^ 
the interfacial free energy. 

It can be shown that the free energy cost of creating a nucleus is related to the probability function P{n), and if the 
appropriate correction is included |28| . the height of the maximum in the curve corresponds to the nucleation barrier, 
while the corresponding cluster size is the size of the critical nucleus. The actual shape of the barrier, as found in 
simulations, will depend on the way in which a cluster is defined. This will therefore also influence the number of 
particles in the critical nucleus. The nucleation barrier itself, however is much less dependent on the cluster definition. 
In the presence of shear the curves have a similar shape. Also here we can identify the top with a critical nucleus size, 
since for smaller clusters there is tendency to shrink, while for larger clusters the growth will dominate. The height 
of the top is a measure for the likeliness that such a cluster will be formed and is therefore related to the time one 
needs to wait in order to observe such a cluster size. 

It is therefore tempting to interpret also the other probability functions of Fig. O in terms of nucleation barriers. 
However, due to the presence of shear in the system, we are dealing with a non-equilibrium situation, and although 
the concept of a nucleation barrier in this case seems intuitively clear, strictly speaking it is wrong since the Gibbs 
free energy is ill defined outside equilibrium. Surprisingly, the application of statistical mechanics to non-equilibrium 
situation can be useful (see e.g, |32|') and it is a canonical challenge to check the validity of carrying over equilibrium 
concepts into non-equilibrium physics. With all the principle caveats which one should have in mind about this, we 
here will assume that we can extend the ideas of classical nucleation theory to the present situation by including 
shear. 



As is the case in classical nucleation theory, we assume that also in the presence of shear a spherical nucleus is 
formed. The two terms in Eq. ^ will therefore remain. It is however reasonable to expect that the difference in 
chemical potential A/x and interfacial free energy jsl will be affected by the shear, but that for moderate shear-rates 
an expansion in powers of the shear rate for both these quantities about their equilibrium values can be performed 



A^ = A^(^«) (1 
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+ 0(7')) 
0(7')) , 



(7) 



where due to the invariance of the shear direction only even powers in the shear rate 7 need to be considered. If we 
insert the expansions ^ in the expression from classical nucleation theory (jSJ we obtain simple expressions for AG* , 
the height of the nucleation barrier and N* , the size of the critical nucleus 
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where the corrections with respect to the unsheared case scale to leading order with 7^. 
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FIG. 4: The height of the nucleation barrier /3AG* as a function of the dimensionless shear rate 'J/k^D for different pressures 
P. The solid lines are parabolic fits through the data. The stars are data from equilibrium Monte Carlo simulations without 
shear. 

In Fig.2|and[Slwe show the results from our simulations where we extracted the height of the nucleation barrier and 
the size of the critical nucleus for various pressures and shear rates. For the zero-shear case, the data from equilibrium 
Monte Carlo simulations have been included. Parabolic fits through the data support the quadratic dependency of 
these quantities on the shear-rate and also enable us to obtain estimates for the second order coefficients in Eqn. [7| 
the results of which are summarised in Table |l| The presence of shear destabilises the solid as is indicated by the 
negative value of cq. 
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TABLE I: Numerical data for different pressures /3P/k'^ on the equilibrium barrier height AG^'^'' , critical nucleus size N^'^''' , 
and second order corrections to the free energy difference and interfacial free energy as obtained from the fitted simulation 
data. 



The shift of the maximum of the curves in Fig. |31as function of the applied shear rate, is interpreted as an increase 
in the nucleation barrier and the critical nucleation size. This indicates that the interaction of the system under 



300 




olOOt 



O 



50 



O pP/K =0.200 
n pP/K^ = 0.224 
O PP/k' = 0.240 



_L 



0.001 



_L 



0.002 



_L 



0.003 



y/K^D 



FIG. 5: The number of particles A^* of the critical nucleus as a function of the dimensionless shear rate ■j/k^D and for different 
pressures P. The solid lines are parabolic fits through the data. The stars are data from equilibrium Monte Carlo simulations 
without shear. 

consideration corresponds to a situation in which crystal nucleation is suppressed by the presence of shear |20j | . In 
addition, the application of shear allows us to supercool the system beyond the point at which in an equilibrium 
simulation spontaneous crystallisation appears. Due to system size restrictions, much larger shear rates could not be 
accessed, because the increase of the critical nucleus size would lead to simulation artefacts, i.e. a percolating cluster. 
It is imaginable that much larger shear rates could enhance crystallisation by pre-layering the system 6o\ . 

It is important to realise that the shear rate should not be considered as a thermodynamic variable 32] . In fact, 
in a recent study of the effect of shear on the location of the solid-liquid coexistence in a Lennard- Jones system, 
Butler and Harrowell found that no purely thermodynamic description of the effect of shear was possible 34] . Shear 
directly affects the transport of particles from the solid to the liquid phase, and this effect is not thermodynamic. The 
expansion in Eq.[7|is simply a way to represent the effect of shear as if it were purely thermodynamic. A priory there 
is no reason to expect the validity of these assumptions and it is therefore very surprising that this simple extension 
of classical nucleation theory is able to describe our observations. 



V. THE STRUCTURE AND RELATIVE ORIENTATION OF THE NUCLEUS 



From an analysis of the bond orientational order parameters it follows that the structure of the nucleus is pre- 
dominantly body-centered-cubic, even though the stable equilibrium phase is the face-centered-cubic crystal. An 
intriguing question is whether the shape and/or orientation of the solid nuclei that are formed in a super-cooled liquid 
are influenced by the presence of moderate shear-rates. Also without shear, a nucleus is in general neither spherical 
nor convex. To this end we calculate for each cluster the moment of inertia tensor. After diagonalisation the principal 
moments of inertia are obtained. For a truly spherical object these should be identical, but due to natural fluctuations 
this is usually not the case. 

In Fig.Elwe show the average ratio A+/Ao and A_/Ao of the largest (A+) and smallest (A_) moment of inertia with 
the intermediate (Aq) one as function of the cluster size n. The dependence on the cluster size indicate that larger 
clusters are slightly more spherical. It is interesting to observe that the shear rates have no visible influence on these 
ratios, which suggest that the shape docs not depend strongly on the applied shear. This is different from the radial 
distribution functions we measured in the liquid under shear, as they become increasingly asymmetric for higher shear 
rates. 

The diagonalisation procedure of the moment of inertia tensor also enable us to characterise the orientation of 
the main axes of the nucleus. Whereas the shape is hardly influenced by the shear, there is a coupling between the 
orientation of the nucleus and the direction of the shear. To quantify this, we calculate the direction of the largest 
principal moment (smallest dimension of the nucleus). The z-component of the corresponding unit- vector gives the 
cosine of the angle with the vorticity direction. This is weighted with the second Legendre polynomial P2 in order to 
take the symmetries into account. The result is shown in Fig.[7| where {P2{z)) is plotted versus the cluster size. For 
an isotropic cluster shape distribution, this correlation function should be zero as is found for the zero-shear case. In 
the presence of shear, however, a negative correlation is found that grows for larger nuclei. This suggests that the 
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FIG. 6: The ratio of both the maximum and minimum moments of inertia with the medium as function of the number of 
particles in the nucleus and for different shear rates and constant pressure j3P/ k? — 0.2. 

smallest dimension of a nucleus have a preference for being oriented perpendicular to the vorticity direction. 
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FIG. 7: The correlation {P2{z)) of the coupling between the direction of the major moment of inertia with the vorticity direction 
z as function of the size of the nucleus and for constant pressure PP/n^ = 0.2. 

In the same manner we can quantify the correlation between the other directions of the principal moments with 
axes of the reference frame. From this analysis it follows that the largest dimension of the nucleus (smallest moment 
of inertia) is preferentially along the vorticity direction, while the smallest dimension prefers to be in the plane of the 
gradient and shear direction. Even within this plane, the distribution of orientations is not isotropic. The measured 
correlations, enable us to calculate the average angle made between the smallest dimension and the gradient direction 
of which the results are shown in Fig.|S| In order to improve the statistics we have averaged over all clusters containing 
more than 100 particles up to the size of the critical nucleus TV*. The inset of Fig. |S1 shows a schematic drawing of 
the preferred orientation of the nucleus with respect to the reference frame and indicates the tilt angle 0. Note that 
the largest dimension of the nucleus (smallest principal moment) is preferably along the vorticity direction, i. e. 
perpendicular to the plane of drawing. 

As can be seen from Fig. |S1 the tilt angle increases linearly with the applied shear rate 7 and only depends weakly 
on the osmotic pressure. Interestingly, a similar tilt occurs when vesicles with a flexible shape are exposed to a linear 
shear flow [3a|. It is, however, unclear whether this is coincidental or not. 



VI. CONCLUSIONS 



In summary we have studied the influence of shear on the nucleation of particles with Yukawa interaction. To this 
end we used a combination of Brownian dynamic simulations with a sort of umbrella sampling technique. The main 
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FIG. 8: The tilt angle 6 of the principal moment of inertia with respect to the y-axis. The inset shows a schematic representation 
of the preferred orientation of the nucleus with respect to the shear direction indicated by the arrows. 

result is that the probability of crystal nucleation decreases with shear rate, while the critical nucleus size increases. 
The shape of nuclei is hardly influenced by the presence of moderate shear rates, but it will affect their orientation. 

In order to understand the qualitative behaviour of the observations, we used an empirical extension of classical 
nucleation theory in which we carry over simple equilibrium concepts to a situation of non-equilibrium. A priory there 
is no reason to expect the validity of such a procedure, but surprisingly it is capable of describing our observations 
and we hope it will stimulate the fundamental research of non-equilibrium physics. 

The observations made here, could in principle be tested experimentally, by studying homogeneous crystal nucleation 
in sheared colloidal suspensions. If a Poiseuille flow as realized in a capillary viscometer 3(| would be used, rather 
than homogeneous Couette flow, we should expect crystal nuclei to appear preferentially in the middle of the flow 
channel. 

For practical reasons we only considered fluids that are weakly sheared, for which one might expect that shear- 
induced layering is unimportant. It is imaginable that for larger shear such pre-ordering phenomena could become 
relevant, and eventually lead to a situation for which crystal nucleation is enhanced |33l |. Another important as- 
sumption has been that we ignored hydrodynamic interactions for computational reasons. Although this is allowed 
for dilute system, this is in general not correct. Closely related is that the implicit solvent is assumed to have a 
homogeneous linear velocity profile and is unaffected by the presence of the more dense crystalline nucleus. For lower 
pressures and higher shear rates this is not correct |37| , as the local transport properties of the liquid can be influenced 
by the nucleus and consequently affect the nucleation process. 

The method can easily be extended to include oscillatory shear [Sg or to consider heterogeneous nucleation [33| near 
a system wall in a sheared suspension. The combination of umbrella sampling with a dynamical simulation method in 
order to obtain information on rare events is generally applicable in equilibrium and non-equilibrium situations and 
is not restricted to Brownian Dynamics. 

Heterogenous nucleation in a sheared system confined between two parallel plates was studied experimentally by 
Palberg and coworkers in another paper of this special issue 19] . A detailed comparison between experiment and 
theory is in principle possible and left for future work. 
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